clc; clear; close all;

% 参数设置
mu = 3.5;
k = 0.25;
y0 = 0;   % 固定 y0 = 0

% 构建初值网格
nx = 500;
nz = nx;
x0_vals = linspace(0, 1, nx);
z0_vals = linspace(-10, 8, nz);

% 迭代参数
N_transient = 2000;  % 瞬态迭代步数
buff_size = 200;    % 用于分类的缓冲区长度
tol = 1e-4;         % 周期检测的容差
max_range = 10;     % 发散检测的范围

% 初始化结果矩阵，basin(i,j) 对应初值 z0_vals(i) 和 x0_vals(j)
basin = zeros(nz, nx);

% 对每个初始条件进行仿真和分类
for i = 1:nz
    for j = 1:nx
        % 设置初始条件 [x0, y0, z0]
        state = [x0_vals(j), y0, z0_vals(i)];
        
        % 瞬态迭代
        for n = 1:N_transient
            [dx, dy, dz] = mclm(state, mu, k);
            state = [dx, dy, dz];
        end
        
        % 利用 buff_size 次迭代构建缓冲区，用于分类 (仅记录 (x,y) 部分)
        buffer = zeros(buff_size, 2);
        for n = 1:buff_size
            [dx, dy, dz] = mclm(state, mu, k);
            state = [dx, dy, dz];
            buffer(n,:) = state(1:2);
        end
        
        % 调用 classify 对当前 (x,y) 序列进行动力学状态分类
        code = classify(buffer, buff_size, buff_size, tol, max_range, mu, k);
        basin(i,j) = code;
    end
end

figure;
imagesc(x0_vals, z0_vals, basin);
set(gca, 'YDir', 'normal');

colormap;

xlabel('x_0');
ylabel('z_0');
title('Basins of Attraction (Enhanced)');